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Abstract 


Flux  limitation  and  preheat  are  important  prooesses  in  eleotron  transport  oocurring  in  laser 
produoed  plasmas.  The  proper  oaloulation  of  both  of  these  has  been  a  subjeet  reeeiving  mueh 
attention  over  the  entire  lifetime  of  the  laser  fusion  projeet.  Where  nonloeal  transport  (instead 
of  simple  single  flux  limit)  has  been  modeled,  it  has  always  been  with  what  we  denote  the 
equivalent  diffusion  solution,  namely  treating  the  transport  as  only  a  diffusion  proeess.  We 
introduee  here  a  new  approaeh  ealled  the  nonloeal  souree  solution  and  show  it  is  numerieally 
viable  for  laser  produeed  plasmas.  It  turns  out  that  the  equivalent  diffusion  solution  generally 
underestimates  preheat.  Furthermore,  the  advanee  of  the  temperature  front,  and  especially  the 
preheat,  ean  be  held  up  by  artifieial  ‘thermal  barriers’.  The  nonloeal  souree  method  of 
solution,  on  the  other  hand  more  aeeurately  deseribes  preheat  and  ean  stably  ealeulate  the 
solution  for  the  temperature  even  if  the  heat  flux  is  up  the  gradient. 
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I  Introduction 


In  a  laser  produced  plasma,  electron  thermal  energy  is  transported  principally  by  electrons 
with  energy  between  about  5  and  16  times  the  thermal  energy.  Electrons  with  less  than  5 
times  the  thermal  energy  do  not  transport  any  signifioant  energy.  Beeause  of  the 
Maxwellian  nature  of  the  electron  distribution  function,  there  are  effectively  no  eleetron  with 
energy  greater  than  16  times  the  thermal  energy.  Thus  even  if  the  mean  free  path  of  thermal 
eleetrons  is  smaller  than  the  temperature  gradient  seale  length,  the  mean  free  path  of  the 
energy  earrying  eleetrons,  (which  is  proportional  to  the  energy  squared)  may  be  longer, 
meaning  the  transport  cannot  be  deseribed  strietly  locally  (for  instance  proportional  to  the 
temperature  gradient). 

This  problem  has  long  been  reeognized  in  laser  produeed  plasmas.  Reeently  we  have 
published  several  papers  [1-3]  examining  the  use  of  a  Krook  model  to  deseribe  nonloeal 
transport.  Of  course  we  are  well  aware  of  the  limits  of  a  Krook  model  and  these  have  been 
fully  diseussed  in  Refs.  [1-3].  There  several  independent  tests  of  the  model  were  made  and 
the  simulations  passed  these  tests.  Also  elements  of  a  fully  eonservative  Krook  model  were 
sketehed  out  [1]. 

Sinee  so  much  of  the  physies  of  laser  matter  interaetion  is  best  modeled  with  fluid 
simulations,  it  is  important  to  find  some  way  of  treating  nonloeal  flux  within  the  constraints 
of  a  fluid  model  if  at  all  possible.  Alternate  simulation  schemes  such  as  Fokker  Planck  or 
direct  simulation  Monte  Carlo  eould,  at  very  great  expense,  treat  nonloeal  transport  more 
aeeurately,  but  would  eertainly  fail  for  treating  every  other  physical  process,  whereas  a  loeal 
fluid  model  works  very  well.  Very  briefly.  Refs.  [1-3]  found  that  the  eleetron  energy  flux  can 
be  split  into  two  parts,  a  diffusive  like  part  describing  the  low  energy  eleetrons,  and  a 
nonlocal  component  deseribing  the  higher  energy  eleetrons.  Regarding  the  loeal  part,  one 
finds  that  the  diffusion  eoefficient  is  redueed  from  elassical.  This  is  the  way  our  Krook  model 
handles  flux  limitation.  The  nonloeal  part  is  expressed  as  an  integral  over  particle  energy  and 
a  eonvolution  over  space.  It  allows  the  hotter  part  of  the  plasma  to  heat  the  cooler  parts  even 
though  these  are  some  distanee  away.  This  is  the  way  our  Krook  model  deseribes  preheat. 

Others  have  also  introduced  variants  of  Krook  models  to  describe  the  nonlocal  transport  [4,5]. 
Previously  other  authors  have  also  proposed  convolutions  to  describe  the  energetie  particle 
transport,  either  ad  hoe  [6],  or  based  on  eomparison  to  Fokker  Planck  solutions  in  infinite 
homogeneous  plasmas  [7,8]. 

What  all  of  these  methods  have  in  common  is  that  when  it  eomes  time  to  perform  a  numerical 
solution  of  the  energy  transport  equation,  they  apply  what  we  term  the  equivalent  diffusion 
solution  [7-9].  That  is  take  the  energy  flux,  however  it  was  derived,  and  divide  by  the 
negative  temperature  gradient  to  get  an  equivalent  diffusion  coefficient  (or  more  aeeurately  an 
equivalent  thermal  eonduetion,  which  has  dimension  nx  /t,  where  n  is  the  electron  density. 
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instead  of  a  diffusion  coefficient  which  has  dimension  x  /t).  Then  the  problem  is  solved  as  a 
diffusion  problem.  This  usually  makes  use  of  an  implicit  numerical  scheme  has  the  advantage 
of  being  unconditionally  stable.  However  in  a  laser  plasma  simulation  the  diffusion 
coefficient  depends  on  temperature  so  the  problem  is  nonlinear  and  more  complicated.  We 
discuss  this  further  shortly.  Furthermore  there  are  problems  with  this  approach  which  were 
apparent  to  its  early  users  [7].  Namely  there  could  be  a  thermal  flux  where  the  temperature 
gradient  is  zero,  implying  infinite  equivalent  diffusion  coefficient.  To  handle  this,  the 
diffusion  coefficient  is  limited  by  some  maximum  value,  typically  a  thousand  or  a  million 
times  classical.  An  even  more  serious  problem  is  that  the  flux  may  be  up  the  temperature 
gradient  rather  than  down  it.  This  would  mean  a  negative  diffusion  coefficient.  The  solution 
(both  analytically  and  numerically)  would  be  unstable.  To  treat  this  case,  wherever  the 
equivalent  diffusion  coefficient  is  negative,  it  is  replaced  with  the  classical  positive  value. 

This  is  the  way  nonlocal  transport  has  been  handled  up  to  now,  but  there  has  neither  been  an 
independent  check  of  it,  nor  examination  of  the  validity  of  the  equivalent  diffusion  solution. 

Here  we  solve  the  nonlocal  transport  equation  without  using  an  equivalent  diffusion  solution 
and  compare  our  solutions  to  it.  We  call  this  new  technique  the  nonlocal  source  solution. 

The  diffusive  part  of  the  transport  is  treated  in  the  usual  way,  with  an  implicit  solution.  The 
nonlocal  part  is  treated  as  an  energy  source  and  treated  explicitly.  One  might  think  there 
could  be  problems  with  numerical  stability.  After  all  the  code  has,  among  other  constraints,  a 
Courant  condition  based  on  the  ion  or  fluid  velocity.  However  as  we  mentioned,  the  electron 
energy  is  transported  by  electrons  with  energy  between  about  5  and  16  times  the  thermal 
energy,  so  one  might  think  such  an  approach  would  be  numerically  unstable.  However  this  is 
usually  not  the  case.  The  implicit  solution  of  the  diffusive  part  gives  rise  to  a  strong 
stabilizing  effect.  In  Section  II  we  work  out  the  numerical  stability  condition  for  the  nonlocal 
source  solution.  Depending  on  the  relative  sizes  of  the  local  versus  nonlocal  flux,  there  may 
be  no  instability  limit  at  all  on  the  time  step.  If  the  nonlocal  part  is  sufficiently  large,  and  the 
local  part  is  sufficiently  small,  there  is  a  stability  condition,  but  the  constraint  is  weak  for 
typical  laser  plasmas. 

It  is  interesting  to  compare  our  Krook  model  for  transport  in  laser  plasmas  with  advanced 
techniques  for  handling  transport  in  magnetically  confined  plasmas  [10,1 1].  For  the  former, 
turbulence  likely  does  not  play  a  major  role,  although  it  may  arise  in  such  a  way  that  it  can  be 
modeled  with  an  enhanced  collision  frequency  [12,13].  Rather  what  is  important  is  the 
nonlocal  transport  of  the  energetic  particles,  the  local  transport  of  the  low  energy  particles, 
and  their  interplay.  Any  theory  here  must  involve  velocity  space.  For  the  magnetic  fusion 
case,  turbulence  apparently  plays  a  major  role,  and  advanced  mathematical  techniques  have 
been  proposed  to  deal  with  it.  However  these  techniques  do  not  involve  velocity  space,  and 
they  are  probably  not  appropriate  for  the  case  of  laser  produced  plasmas. 

Section  III  works  out  the  stability  condition  for  our  Krook  model  of  electron  transport.  This 
lets  us  evaluate  the  maximum  stable  time  step  and  insure  that  the  code  runs  with  a  smaller 
time  step.  In  most  cases  we  have  looked  at,  there  is,  in  practice,  no  additional  constraint. 
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Section  IV  applies  this  to  a  simple  problem,  a  uniform  (i.e.  infinitely  massive  ions)  plasma 
with  a  background  temperature  of  1  keV,  but  imbedded  within  is  a  slug  of  50  keV  plasma. 

We  solve  for  the  time  evolution  of  the  temperature  profile  using  classical  thermal 
conductivity,  Krook  model  using  an  equivalent  diffusion  solution,  and  Krook  using  a  nonlocal 
source  solution.  For  the  latter  we  find  that  the  numerical  stability  condition  is  well  satisfied. 
The  classical  solution  has  a  rapidly  moving  heat  front  and  no  evidence  of  preheat,  as 
expected.  The  equivalent  diffusion  solution  has  a  slower  moving  main  heat  front  and  a  small 
amount  of  preheating.  The  nonlocal  source  solution  has  a  slower  moving  heat  front  than 
classical,  but  faster  than  the  equivalent  diffusion  solution.  Most  important,  it  has  stronger 
preheat  than  the  equivalent  diffusion  solution.  Our  conclusion  is  that  the  equivalent  diffusion 
solution  can  underestimate  the  preheat.  Also  we  show  that  the  nonlocal  source  solution  has 
no  problem  with  energy  flux  up  a  temperature  gradient,  whereas  the  equivalent  diffusion 
solution  encounters  artificial  ‘thermal  barriers’.  Section  V  briefly  reexamines  the  problem  of 
a  4  keV  plasma  slug  embedded  in  a  1  keV  plasma.  This  was  initially  examined  with  a  Fokker 
Planck  simulation  by  Matte  and  Virmont  [14].  In  Ref.  (2)  we  compared  the  equivalent 
diffusion  solution  to  this  problem  with  the  Fokker  Planck  solution,  and  concluded  that  the 
Krook  model  compared  more  favorably  than  any  other  transport  model,  especially  in  the  time 
asymptotic  limit.  In  Section  V.  we  examine  this  problem  comparing  the  nonlocal  source 
solution  with  the  equivalent  diffusion  solution  as  a  function  of  time.  We  find  the  two 
solutions  are  close,  and  in  the  time  asymptotic  limit  approach  one  another.  Finally  Section  VI 
gives  conclusions. 


II.  The  nonlocal  source  term  approach  and  its  stability 

In  the  models  for  non  thermal  electron  energy  transport  which  we  have  developed  [1-3],  we 
found  that  the  electron  thermal  energy  flux  is  the  sum  of  two  terms,  a  local  term  which  has  the 
properties  of  classical  transport,  and  a  nonlocal  term  which  is  a  convolution  over  space  as 
well  as  an  integral  over  energy.  As  the  temperature  equation  has  many  terms,  representing 
many  physical  processes  in  it,  we  examine  the  numerical  stability  of  a  simpler  system,  namely 
a  temperature  equation  of  the  form 


dt 


=  WD 


d^T 


+ 


j  Js||-jjV^«(s)exp{-|A(S)(x- 


(1) 

ox 


For  an  actual  laser  plasma,  the  situation  is  more  complicated  because  D  is  a  function  of  T(x) 
and  a  is  a  function  of  T(x’).  However  we  proceed  with  the  simple  linear  case  and  discuss  its 
validity  and  applicability  shortly.  In  an  analytic  theory  of  course,  the  nonlocal  flux  alone  is 
sufficient.  But  when  interpolating  onto  a  finite  spatial  grid,  for  many  energies  (i.e  S’s)  the  x’ 
integration  takes  place  within  a  single  grid  cell.  For  these  energies,  the  finite  difference 
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approximation  is  not  even  qualitatively  correet.  Thus  in  practice,  it  is  necessary  to  split  the 
flux  into  a  local  and  nonlocal  part.  We  will  discuss  this  in  more  detail  shortly. 

Here  D  is  the  classical  diffusion  coefficient,  and  W  is  the  reduction  from  the  classical  result 
[1-3].  It  arises  because  as  nonlocal  effects  become  more  important,  the  local  transport 
decreases.  The  second  term  in  Eq.  (1)  comes  from  the  thermal  flux  arising  from  nonlocal 
effects.  It  is  an  integral  over  particle  energy  S.  As  discussed  in  [2,3],  the  minimum  energy 
Scr  is  determined  first  by  the  condition  that  WD  is  positive  (meaning  Scr/T>3.5);  and  second 
by  the  condition  that  the  integral  over  x’  is  reasonably  represented  by  the  summation 
(meaning  K(Scr)Ax<0.2;  K(S)  being  a  decreasing  function  of  S  reflecting  the  longer  mean 
free  paths  of  the  more  energetic  particles).  The  last  terms  on  the  right  simply  defines  this  rate 
of  change  of  temperature  in  terms  of  the  gradients  of  local  and  non  local  fluxes.  As  long  as 
WD  and  a  are  positive,  there  are  no  unstable,  i.e.  exponentially  growing  solutions  to  Eq.  (1). 

In  the  literature  [1-9],  this  equation  has  typically  been  solved  with  what  we  call  the  equivalent 
diffusion  solution.  Namely  one  defines  an  equivalent  diffusion  coefficient 


D 


eq 


(ll  +  (Inl 

dT 

dx 


(2) 


and  solves  the  equation  as  a  diffusion  equation.  However  if  the  temperature  gradient  becomes 
zero,  Deq  is  infinite.  Even  more  of  a  problem,  if  Dgq  is  negative  in  some  regions,  the  diffusion 
equation  is  violently  unstable.  To  treat  these  cases,  the  approach  has  been  to  restrict  Deq  to 
some  maximum,  say  10^  times  classical,  and  in  regions  where  Dgq  is  negative,  to  replace  it 
with  the  positive  classical  result.  We  will  see  that  this  can  lead  to  greatly  inaccurate  solutions 
in  some  cases. 

Instead  of  the  equivalent  diffusion  solution,  we  treat  the  added  term  in  Eq.  (1)  as  a  nonlocal 
source  term  and  treat  it  explicitly,  that  is  in  solving  for  the  (n-Hl)  time  step,  the  source  term  is 
evaluated  at  the  n^’’  time  step.  However  we  treat  the  local  diffusion  term  implicitly,  that  is  the 
diffusion  term  is  evaluated  at  the  (n  -I-  1)  time  step.  If  a  diffusion  equation  is  solved 
implicitly,  there  is  no  constraint  on  the  time  step  for  numerical  stability,  but  only  for 
accuracy.  However  when  treating  diffusion-like  equations  explicitly,  the  condition  for 
numerical  stability  demands  a  very  small  time  step.  Indeed  for  the  case  of  a  grid  of  space  step 
Ax  and  time  step  At,  the  well  known  condition  for  stability  is  (for  W=l),  At<Ax  /2D. 

Let  us  pause  here  to  briefly  discuss  the  implicit  solution  to  the  temperature  equation  as  done 
in  laser  fusion  calculations.  If  D  were  constant,  the  implicit  calculation  is  straightforward. 

The  problem  is  that  D  has  strong  temperature  dependence;  it  varies  roughly  as  T  ,  with  other 
weaker  dependences  as  well.  Each  laboratory  doing  laser  fusion  calculations  uses  its  own 
approach.  At  NRL,  we  use  the  fluid  simulation  code  FAST  [15].  It  treats  the  implicit 
temperature  calculation  with  an  iteration  scheme  [16].  The  temperature  dependence  of  D  is 
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expressed  at  the  current  time  step  and  then  T  is  advanced  implicitly  as  a  linear  calculation. 
This  new  T  is  then  used  in  D  and  the  process  is  iterated  until  convergence  is  achieved.  In 
practice,  it  works  very  well  with  a  few  iterations. 

Hence  before  applying  a  nonlocal  source  term  approach,  we  must  first  examine  the  numerical 
stability  where  the  source  term  is  treated  explicitly  and  the  diffusion  term  is  treated  implicitly. 
An  actual  numerical  stability  calculation,  taking  into  account  the  dependence  of  D  and  a  on  T 
would  be  extremely  complicated  due  to  the  fact  that  this  introduces  both  inhomogeneity  and 
nonlinearity.  Furthermore,  even  if  one  could  do  such  a  stability  calculation,  it  would  be  of 
little  use;  every  temperature  profile  would  have  a  different  stability  condition.  To  proceed, 
we  do  a  much  simpler  numerical  stability  calculation,  we  consider  D  and  a  to  be  independent 
of  space  (but  a  depends  on  S).  Hence  we  consider  only  the  linear  homogeneous  problem. 
Intuitively  one  would  expect  this  to  give  reasonable  guidance;  such  an  approach  is  standard  in 
the  literature  and  goes  back  to  a  textbook  on  numerical  methods  [17].  Furthermore,  we  will 
see  shortly  that  there  are  other  convincing  reasons  to  believe  that  the  results  of  the  stability 
analysis  can  be  applied  in  a  much  more  general  sense. 

To  continue,  we  write  out  the  discrete  forms  of  the  various  operators.  Here  superscripts  refer 
to  time  step,  subscripts  to  spatial  step.  To  examine  the  numerical  stability,  we  assume  that 

=  [exp(- v'At)]r"p  (i.e.  T(t)  ~  exp(-vt))  (3a) 


r"p+i  =  [exp(-A:Ax)]r"p  (i.e.  T(x)~exp(ikx))  (3b) 

Ultimately  we  wish  to  solve  for  the  dispersion  relation  between  v  and  k.  Note  that  if  v<0,  for 
the  finite  difference  problem,  whereas  if  v>0  for  the  continuum  problem  (as  it  is  for  our  case), 
the  solution  is  numerically  unstable. 

We  now  write  out  the  discrete  forms  of  the  various  terms  in  Eq.  (1).  The  first  is 

dt  At 

The  second  is 


exp(ikAx) -2  +  exp(-ikAx) 

dx  1 


T"p  exp(-vA^) 


(5) 


where  we  note  here  that  on  the  right  hand  side  of  Eq.  (5),  T  is  evaluated  at  the  advanced  time 
step,  t=(n+l)At.  The  nonlocal  term  on  the  right  hand  side  of  Eq.  (1),  when  written  in  finite 
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difference  form,  and  using  Eq.  (3b)  breaks  up  into  summations  over  various  geometric  series 
with  ratio  exp[(-K(S)+ik)Ax]  or  exp[(-K(5)-ik)Ax].  Doing  these  sums  over  the  geometric 
series,  we  find  that  this  term  is 

-  =  1"  p\d’R  ^^^^{exp(ikAx)  -  2  +  exp(-/kAx)} 

dx  •'  2Ax 

^  1  ^  exp[(-if(5)-iA:)Ax] 

1  -  exp(-  ^(5) Ax  +  /kAx)  1  -  exp(-  if  (5) Ax  -  ikAx) 

=  T" p  [  ^^""^‘^^""^{exp(/A:Ax)  -  2  +  exp(-  iA:Ax)}i?(^(S)Ax,A:Ax))  (6) 

•'  2Ax 

Hence  we  find  that  the  dispersion  relation  is 

=  {exp(,mv)-  2  +  exp(-  ,*Ax)[ +  f<i5  ;;(g(S)A,v.tAx)) 

(7) 

As  is  usually  the  case  in  investigations  of  numerical  stability  [17],  and  as  we  find  by  an 
investigation  of  Eq.(7),  the  most  unstable  wave  number  satisfies  kAx=7T:.  To  show  this,  the 
only  term  in  Eq.  (7)  which  can  drive  instability  is  the  R  term,  the  nonlocal  term.  The 
numerical  solution  will  first  exhibit  instability  where  this  term  maximizes.  By  manipulating 
the  expression  for  R,  and  using  the  half  angle  formula,  it  is  possible  to  show  that  each  term  in 
the  S  summation  has  for  its  kAx  dependence 

1  -  cos  kAx 
x4(S)  -  cos  kAx 

where  A(S)  is  a  function  of  S  which  is  not  necessary  to  write  out  for  our  purposes  here 
except  to  note  that  is  always  greater  than  one.  Clear  this  maximizes  at  kAx  =ti,  meaning  that 
the  numerical  scheme  is  first  unstable  for  this  wave  number. 

Eet  us  now  reconsider  the  linear  homogeneous  approximation  used  in  the  stability  analysis. 
Even  though  the  temperature  varies  in  x,  there  is  rarely  very  much  temperature  variation 
within  a  few  grid  cells.  However  numerical  instability  occurs  on  the  very  shortest  spatial 
scale  length,  so  the  background  appears  nearly  homogeneous.  Thus  the  approximation  of 
constant  D,  where  the  local  temperature  used  in  D  is  regarded  as  a  constant  local  value  is 
reasonably  valid. 

Taking  kAx=7T:,  we  find 
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(8) 


i?(^(S)Ax:,;r) 


1  -  exp(- ^(S)Ax:) 
1  +  exp(- ir(S)Ax) 


KiE)Ax 

2 


We  assume  here  that  k(S)Ax  «1.  For  those  energies  for  which  kAx»l,  expressing  the  dx’ 
integral  in  Eq.  (1)  as  a  finite  summation,  as  in  deriving  Eq.(6)  is  not  even  qualitatively 
correct.  For  these  energies  for  which  k(S)Ax  »1,  the  energy  flux  is  treated  with  the  local 
approximation.  In  treating  the  energy  flux  in  for  instance  a  laser  produced  plasma,  there  is  a 
very  wide  range  of  particle  energies,  say  from  S=  lOeV  to  S  =  lOOkeV.  Since  the  mean  free 
path,  proportional  to  S  '  varies  over  8  orders  of  magnitude,  and  since  a  laser  plasma 
simulation  typically  has  fewer  than  1000  spatial  cells,  flux  cannot  be  treated  either  all  locally 
or  all  nonlocally  (unless  the  most  energetic  electron  has  mean  free  path  less  than  a  grid  cell). 

Let  us  illustrate  this  point  with  a  simple  example.  Say  the  plasma  temperature  is  1  keV  and 
the  mean  free  path  of  a  2  keV  electron  is  one  grid  cell.  In  practice,  we  find  that  due  to  the 
Maxwellian  nature  of  the  distribution  function,  there  is  no  contribution  to  the  thermal  flux 
from  energies  above  16  times  the  temperature  [1].  Since  the  mean  free  path  is  proportional  to 
the  temperature  squared,  this  means  that  the  most  energetic  electron  spreads  out  its  energy 
over  64  grid  cells.  Clearly  this  must  be  treated  nonlocally.  On  the  other  hand,  there  are  also 
many  electrons  with  energy  of  for  instance  500  eV.  These  deposit  their  energy  in  1/16  of  a 
grid  cell.  Clearly  any  numerical  algorithm  which  has  them  spreading  their  energy  out  over 
many  grid  cells  will  be  incorrect.  Thus  given  the  constraints  of  a  finite  grid,  one  has  no 
choice  but  to  break  up  the  energy  range  into  two  regimes,  the  lower  energies  treated  locally, 
the  higher,  nonlocally. 

We  find  that  a  good  approximation  is  to  treat  the  electron  flux  nonlocally  if  its  mean  free 
path  is  five  grid  cells  (i.e  kA^<0.2)]  or  longer  [3]. 

So  the  dispersion  relation  for  this  kAx  is 


exp(-  vAt)  = 


l-AAt 


\  +  AWDAt/ Ax 


2  ’ 


with 


^  =  I  dRK^  (S)a(S) 


(9) 


The  condition  for  numerical  instability  is  that  the  right  hand  side  of  Eq(9)  has  a  magnitude 
greater  than  unity.  Since  the  sign  in  this  case  is  negative,  this  means  that  for  numerical 
instability,  we  need  Re  v  <0  and  Im  v  =7r.  Hence  the  condition  for  numerical  instability  is 


AAt  >2  + 


AWDAt 

Ax" 


(10) 


2 

As  long  as  A<4WD/Ax  ,  solving  the  local  diffusion  part  implicitly,  and  the  nonlocal  source 
term  explicitly,  will  always  be  numerically  stable.  In  other  words  the  natural  stability  of  the 


implicit  solution  of  the  solution  of  the  diffusion  equation  will  always  overpower  the  tendeney 
of  the  explieit  solution  of  the  nonloeal  souree  term  toward  stability.  This  is  a  numerieal 
condition  involving  only  Ax,  and  not  At.  On  the  other  hand,  if  A>4WD/Ax  ,  we  find 
numerieal  stahility  as  long  as  the  time  step  satisfies  the  condition 


^  1 

AM^J2-W 


(11) 


2 

Here  Atgx  is  the  maximum  time  step  for  stahility  for  the  explicit  diffusion  calculation.  Ax  /2D. 
This  is  plotted  graphieally  in  Fig  (1).  In  this  case,  the  stability  condition  once  again  depends 
on  At. 


III.  Caleulation  of  a(S)  and  the  numerieal  stahility  eondition  for  a  laser  produced  plasma. 

In  our  earlier  work,  we  have  written  out  expressions  for  the  thermal  flux  using  a  Krook  model 
for  eleetron  eleetron  collisions  and  a  Fokker  Planck  model  for  eleetron  ion  collisions  [1-3]. 
There  the  eleetron  ion  and  electron  collision  frequencies  were 


= 


5.8x10  «A 


2v„ 


,  3/2 


1  +  1 


(5.8/7.7)(s/r  f ']  [l  +  (5.8/7.7)(S/r)''" 


l^v,SiE,Z)  (12) 


v^,=L34ZvXmv^ ^v^S,iE,Z) 


(13) 


Vei  +  Vee  =  VeS(v,Z)  (14) 

Notice  that  these  eollision  frequencies  are  funetions  of  particle  energy  S,  as  well  as  position 
through  the  spatial  dependenee  of  the  eleetron  density  n,  the  Coulomb  logarithm  A,  the 
eleetron  temperature  T,  and  the  charge  state  Z. 


In  Refs.  [1-3]  we  derived  an  expression  for  the  energy  flux: 


qix)  =  -KX'"T 


sp 


/lioo) 


Z^jdiyfA) 


(x) 

2;^(oo) 


[I] 

[l  0.92 

J  J 

-cc 

Z2(x)S 

(x,x',S))  =  (x)  +  (x) 


(15) 
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where  qL  and  Qnl  are  the  local  and  nonlocal  terms  in  q(x)  and 


.A(x',5)S^ 

S  Z4(x) 

[exp{-  [s/r  (x')]}exp{-  i/(x,x'5)} 

T  (x')  Z2(x') 

I  dUx') 

A(r’)r^(r’)''"{l  +  Z(r')/2} 


(16) 


Other  quantities  in  Eq.  (15)  are 


p{y)  = 


\u  /  Z2 } 
S{u) 


exp(-w) 


R{y)  = 


S{u) 


exp(-w) 


(17) 


(18) 


Z.(Z)  =  juW5^^h  (19) 

i  S(u,Z) 

Analytic  approximations  to  these  quantities  have  been  given  elsewhere  (1-3).  We  note  that  in 
the  analytic  approximations,  the  Z  dependence  is  only  meaningful  if  Z>1 .  However  the 
colder  regions  of  a  laser  produced  plasma  are  not  necessarily  fully  ionized.  In  these  regions, 
we  use  the  analytic  expressions,  but  evaluate  them  at  Z=l.  Regarding  the  nonlocal  transport, 
these  cold  unionized  or  partially  ionized  regions  play  no  role  as  a  source  of  energetic 
electrons,  but  of  course  they  can  be  important  as  regards  other  physical  effects. 

The  quantity  H  is  given  by 


H{x,x\K)  =  ^dx"K{x'\K) 

x' 


where 


V {r,  S)[v^^  (r,  S)  +  (r ,  S) J 

K{r,id)  = - 


6x10 'S 


1/2 


And  Ksp  is  the  Spitzer  conductivity 


(20) 


(21) 
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(22) 


K 


sp 


3.1x10'°  0.4(Z  + 0,235) 
A  Z  +  4.242 


and  the  expression  for  from  Refs.  (2  and  3)  is 


(x)  =  Max<  5T^  {x),Min 


(1 . 1 5x1  O'"  n(x)A(x)(l  +  Z(x)  /  2)Ax)""  , 
(2.03x10 n(x)A(x)(l  +  Z(x)  /  2% 


(23) 


In  Equations  (12  -23),  all  quantities  are  in  egs  units  except  energies,  (5  and  T)  which  are  in 
electron  volts.  Since  the  flux  is  the  energy  flux,  to  write  the  equation  in  terms  of  diffusion 
equation  for  temperature,  we  must  divide  the  flux  by  the  specific  heat  Cv,  which  in  the  case  of 
a  fully  ionized  plasma  is  3n/2. 

Hence  the  dimensionless  quantity  W  is  given  by 


W(EM/TXx))  = 


As,.(x)/r(x)) 


^  ^^^Z,R{EMIT{x)) 


(24) 


As  shown  in  Ref.  (3),  W(u)  need  be  evaluated  only  for  5  <  u<  14.  It  varies  (monotonically 
increasing)  from  about  0.125  at  u=5  to  unity  for  u=14  and  above.  Because  W<1  where 
nonlocal  transport  occurs,  the  flux  due  to  the  ‘local’  transport  is  reduced. 

As  explained  in  Refs.  (1-  3),  Scr  is  the  critical  energy  dividing  the  low  energy  electrons,  with 
short  mean  free  path,  which  must  be  treated  locally;  from  the  higher  energy,  long  mean  free 
path  electrons  which  must  be  treated  nonlocally.  As  discussed  in  Refs.  (2  and  3),  this  depends 
on  the  size  of  the  grid  cell  Ax.  This  is  a  consequence  of  the  rather  small  number  of  total  grid 
cells  (several  hundred)  to  describe  a  range  of  mean  free  paths  extending  over  a  range  of  a 

o 

factor  of  more  than  10  in  lengths.  However  there  is  an  additional  consideration,  namely 
Scr/T  >  5.  We  discuss  further  the  consequences  of  this  choice  in  the  Appendix. 

A  laser  produced  plasma  is  extremely  inhomogeneous.  However  in  the  stability  calculation  in 
the  last  section,  we  assumed  a  homogeneous  system.  But  the  numerical  stability  is 
determined  at  the  shortest  wavelength,  kAx=7r.  Hence  as  regards  the  stability  calculation,  we 
may  consider  the  system  to  be  homogeneous,  so 


H=K(x,5)[x-x’] 


(25) 


and 
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(26) 


r 

Y 

exp 

2 

7  T 

1  0  9"’'^"^ 

N 

ts) 

1  _ 

f  T) 

3« 

/3(cc) 

1  \)  ,y  z. 

7  ” 

L  -^2^  J 

^5/2|-i  + 

Z/2] 

V 

A 

Here  all  quantities  are  defined  as  a  function  of  x  (recall  Scr  is  a  function  of  x)  and  5. 

To  find  A,  the  quantity  used  in  the  stability  condition,  the  integral  of  K  (S)a(S)  must  be 
carried  out  from  Scr  to  infinity.  Since  is  proportional  to  S'"^,  this  integral  is  particularly 
simple  to  perform  analytically;  however  to  get  a  simple  analytic  expression,  one  must  make 
an  asymptotic  approximation  to  the  part  of  the  total  integral  which  goes  as  x'^e'^'. 


Doing  the  S  integrals,  we  find 


A  = 


p{^)T 


5/2 


-expl 


( 

S  ^ 

— 

cr 

< 

V 

T  J 

^  +  1-1.9^  +  0.921 

T  Z, 


( Z 

'^4 

2 

T 

7i ) 

M 

^  cr 

(27) 


Thus  using  the  spatially  dependent  expressions  for  A  and  W  appropriate  for  a  laser  produced 
plasma,  we  can  determine  the  maximum  stable  time  step  for  stability  in  any  such  simulation. 
This  is  then  added  to  the  constraints  on  the  time  step. 


IV.  A  sample  calculation:  Failure  modes  of  the  equivalent  diffusion  solution 

Here  we  compare  the  nonlocal  source  term  solution  and  the  equivalent  diffusion  solution  for  a 
simple  test  problem.  A  uniform  deuterium  plasma,  with  electron  density  8x10  cm'  andZ=l 
is  set  up  over  a  region  of  0.5  cm.  On  the  leftmost  0.05  cm,  the  temperature  is  55  keV.  Over 
0.03  cm,  the  temperature  drops  to  IkeV.  For  x>0.08  cm,  the  initial  temperature  is  constant  at 
IkeV.  A  increases  monotonically  with  temperature  from  3.4  at  1  keV  to  7.3  at  55  keV.  This 
temperature  profile  is  set  up  at  t=0.  The  left  boundary  is  maintained  at  55  keV,  while  the 
temperature  gradient  is  taken  as  zero  on  the  right  boundary.  If  the  nonlocal  expression  for  the 
flux  here  is  not  zero,  the  equivalent  diffusion  model  would  give  an  infinite  diffusivity. 
However  this  has  little  effect  on  the  calculation  and  is  treated  by  the  maximum  D  discussed  in 
the  introduction.  A  much  more  complete  discussion  of  boundary  conditions  is  given  in  Ref. 
[3].  The  time  and  space  dependence  of  the  temperature  is  followed  using  a  variety  of 
transport  models.  This  problem  is  rather  like  that  examined  by  Matte  and  Virmont  [14],  and 
discussed  in  [2],  except  that  here,  there  is  a  much  greater  range  of  temperature  variation,  i.e. 
over  a  factor  of  50.  The  mean  free  path  (i.e.  K'^)  of  a  100  keV  electron  is  about  0.01 1  cm. 

The  temperature  diffusion  equation  is  solved  by  breaking  the  region  into  a  grid  of  128  spatial 
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steps.  In  all  cases,  we  find  that  the  time  step  is  stable  as  specified  in  Section  II  using  the 
evaluations  of  A  and  W  as  spelled  out  in  Section  III. 

The  solution  using  pure  Spitzer  (classical  conductivity)  is  shown  in  Fig  (2).  Because  the 
thermal  diffusion  coefficient  is  proportional  to  the  implicit  procedure  is  not  as  simple  as 
indicated  in  Section  II.  This  temperature  must  also  be  evaluated  at  the  advanced  time  step, 
and  there  is  no  simple  way  to  do  this.  Accordingly  all  implicit  schemes  in  laser  plasma 
simulations  initially  take  T  at  the  previous  time  step  where  it  is  known,  and  use  an  iteration 
process  to  converge  to  the  proper  implicit  result.  The  time  interval  between  the  different 
temperature  profiles  is  600  ps.  In  Fig  (3a)  is  shown  the  corresponding  plots  for  the  solution 
using  the  equivalent  diffusion  solution.  Because  W  is  less  than  unity,  most  of  the  front  moves 
slower  than  it  does  for  the  classical  case.  However  there  is  some  preheating  that  is  apparent. 

In  Fig  (3b)  is  shown  the  solution  using  the  nonlocal  source  solution.  Again  the  main  front 
moves  outward  noticeably  slower  than  in  the  classical  calculation,  but  faster  than  for  the 
equivalent  diffusion  solution.  However  now  there  is  noticeably  more  preheat.  One  might 
think  that  in  the  equivalent  diffusion  solution,  we  simply  did  not  let  the  maximum  diffusion 
coefficient  be  large  enough,  but  this  is  not  the  case.  We  did  this  calculation  where  the 
maximum  diffusion  coefficient  was  limited  to  10  ,  10  ,  and  10  times  classical,  and  there  was 
virtually  no  difference  between  the  results. 

The  reason  for  so  much  less  preheat  in  the  equivalent  diffusion  calculation  seems  to  be  that 
here,  the  preheating  is  modeled  as  a  diffusion  process,  so  to  get  from  point  a  to  point  b,  the 
energy  has  to  pass  through  all  points  in  between  at  whatever  the  maximum  diffusion  speed  is. 
However  in  the  nonlocal  source  solution,  heating  from  the  energetic  electrons  is  virtually 
instantaneous,  at  least  on  the  time  scale  of  the  electron  electron  collision  time  (~10'  sec  for  a 
100  keV  electron).  Thus  the  equivalent  diffusion  solution  does  deemphasize  preheat  as 
compared  to  the  nonlocal  source  solution. 

The  role  of  preheat  can  be  illustrated  more  clearly  by  examining  a  third  possible  solution.  In 
Fig.  (3  c)  is  shown  the  result  of  the  nonlocal  source  solution,  but  where  the  nonlocal  part  is 
itself  artificially  set  to  zero.  This  is  equal  to  the  equivalent  diffusion  solution  where  the 
nonlocal  source  is  also  set  equal  to  zero.  We  see  that  as  in  the  Spitzer  case,  there  is  no 
preheat  at  all,  and  the  main  heat  front  advances  even  more  slowly  than  in  Fig.  (3a).  Thus  the 
following  rather  complicated  picture  emerges.  In  all  cases,  the  main  front  advances  more 
slowly  than  it  does  for  classical  thermal  conduction.  This  is  due  to  flux  limitation, 
characterized  by  the  fact  that  W(u),  from  Eq.  (23)  is  less  than  unity.  We  can  then  regard 
going  from  Fig  (3c  to  3a  to  3b)  as  adding  more  and  more  preheat.  We  see  that  the  effect  of 
preheat  is  twofold.  First  there  is  the  preheat  itself  However  secondly,  the  preheat  increases 
the  temperature  just  ahead  of  the  temperature  front.  This  increases  the  ‘local’  flux  there  and 
speeds  up  the  advance  of  the  main  front.  Thus  preheat  contributes  significantly  to  pulling  the 
main  temperature  front  forward. 
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Another  failure  mode  of  the  equivalent  diffusion  model  ean  be  shown  even  more  strikingly 
and  clearly.  It  results  from  the  fact  that  where  the  flux  is  up  the  gradient,  the  equivalent 
diffusion  model  replaces  the  flux  with  the  classical  flux,  which  may  be  much  smaller  and  has 
the  opposite  sign.  Instead  of  an  initially  uniform  temperature  of  IkeV  from  0.08  to  0.5  cm,  let 
the  temperature  decrease  linearly  in  space  from  IkeV  at  0.08  cm  to  500  eV  at  0.3  cm,  and 
then  rise  linearly  to  IkeV  again  at  0.5  cm.  As  far  as  the  nonlocal  flux  is  concerned,  it  is 
outward,  and  barely  notices  the  small  temperature  dip  and  rise.  However  the  equivalent 
diffusion  solution  sees  a  flux  up  the  temperature  gradient  for  x>0.3  cm,  and  to  maintain 
numerical  stability,  it  replaces  the  nonlocal  flux  up  the  temperature  gradient  with  the  much, 
much  smaller  classical  flux  down  the  gradient.  In  the  equivalent  diffusion  model  this 
temperature  dip  presents  a  thermal  barrier,  the  heat  cannot  get  to  the  other  side  without  first 
diffusing  through  a  region  where  the  temperature  gradient  is  very  small,  i.e.  before  filling  in 
the  dip. 

Figure  (4a  and  b)  show  the  temperature  profiles  as  a  function  of  time  for  the  equivalent 
diffusion  solution  (4a)  and  for  the  nonlocal  source  solution  (4b).  Notice  that  the  initial  dip 
and  rise  in  temperature  is  hardly  visible  on  the  scale  shown.  However  the  equivalent 
diffusion  solution  hangs  up  around  the  temperature  minimum  whereas  the  nonlocal  source 
solution  does  not.  This  is  even  more  apparent  in  a  blown  up  temperature  scale  shown  in  Fig. 
(5a  and  b).  The  equivalent  diffusion  solution  (5a)  clearly  hangs  up  around  the  temperature 
minimum  and  cannot  get  by  it  for  quite  some  time.  The  nonlocal  source  solution  (5b) 
however  easily  overrides  the  temperature  dip  and  gets  past  it  in  a  very  short  time. 

Notice  that  even  though  the  energy  flux  is  going  up,  rather  than  down  a  temperature  gradient, 
the  calculation  is  stable.  This  is  as  determined  in  Section  II;  we  have  checked  that  the  time 
step  is  always  below  the  stability  threshold. 

To  see  in  more  detail  the  effect  of  this  thermal  barrier,  we  compare  the  equivalent  diffusion 
solution  for  the  case  where  there  is,  and  is  not,  a  temperature  dip.  Figure  6  shows  blow  ups 
plots  of  the  temperature  profile  at  times  0.29,  1.49,  2.69,  and  3.88  ns.  The  dotted  curves  are 
calculations  with  no  initial  temperature  dip;  the  solid  curves  are  as  in  Fig.  5.  The  main 
temperature  front  advances  somewhat  more  slowly  for  the  case  with  the  initial  temperature 
dip.  However  the  greatest  difference  is  in  the  preheat.  For  no  initial  temperature  dip,  we  see 
significant  preheat,  but  not  nearly  as  much  as  for  the  equivalent  source  solution.  However  for 
the  case  of  the  initial  temperature  dip,  where  on  the  right  hand  side  of  the  dip,  the  diffusion  is 
classical  and  to  the  left,  there  is  virtually  no  preheat. 

To  test  the  accuracy  of  our  numerical  scheme,  we  reran  the  case  of  Fig  (3b),  the  nonlocal 
source  solution,  on  two  different  spatial  grids.  First  we  took  our  standard  case  of  128  grid 
cells,  and  then  we  took  256  grid  cells.  Here  the  maximum  temperature  was  taken  as  45  keV 
instead  of  55.  A  plot  of  the  position  of  the  advancing  front,  defined  as  the  position  of 
maximum  curvature  of  the  T(x)  graph,  as  a  function  of  time  is  shown  in  Fig  7.  Clearly  to 
within  an  accuracy  of  a  few  percent,  the  results  are  independent  of  grid  size. 
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To  summarize,  the  equivalent  diffusion  solution  unphysically  reduces  the  effect  of  preheat  as 
compared  to  the  more  accurate  nonlocal  source  solution.  Furthermore,  it  is  susceptible  to 
additional  inaccurate  results  because  artificial  thermal  barriers  can  be  set  up. 


V.  A  reexamination  of  the  Matte  Virmont  Fokker  Planck  solution. 

In  Ref  (2)  we  compared  the  Krook  solution  to  that  of  Matte  and  Virmont  [14].  They  used  a 
Fokker  Planck  code  to  examine  the  evolution  of  a  temperature  profile  which  starts  out  at  4 
keV  on  the  left  boundary  and  drops  sharply  to  1  keV.  Initially  the  1  keV  temperature  is 
constant  to  the  right  hand  boundary.  For  the  plasma  in  Ref.  (10),  Z=6,  the  electron  density  is 
n=  6x10  cm'  ,  the  time  t  is  measured  in  collision  times  for  a  particle  at  the  average 
temperature,  and  the  length  L  is  measured  in  units  of  mean  free  path  lengths  for  this  particle. 
The  particular  simulation  done  here  is  for  an  L=20,  which  for  our  system  is  0.00165  cm. 
Figure  8  shows  the  results  of  the  equivalent  diffusion  solution  (solid)  and  nonlocal  source 
solution  (dashed)  as  the  time  increases.  The  single  solid  curve  is  the  initial  temperature. 
Clearly  the  two  solution  methods  are  reasonably  close  to  each  other  for  this  problem,  but  in 
the  transient  regime,  the  nonlocal  source  solution  shows  somewhat  more  preheat.  In  the  time 
asymptotic  limit,  the  two  solutions  are  virtually  the  same  and  match  well  to  the  time 
asymptotic  Fokker  Planck  solution  as  discussed  in  Ref.  (2). 


VI.  Conclusions 

In  using  a  nonlocal  transport  model,  we  find  that  it  is  not  necessary  to  use  an  equivalent 
diffusion  solution  to  the  energy  transport  equation.  The  new  method,  called  the  nonlocal 
source  solution,  relies  on  splitting  the  expression  for  the  heat  flux  into  a  Spitzer  like  term, 
proportional  to  minus  the  temperature  gradient,  and  a  nonlocal  source  term.  We  have 
investigated  the  numerical  stability  and  found  that  it  usually  allows  time  steps  comparable  to 
time  step  constraints  already  in  the  code.  Thus  using  a  nonlocal  source  solution,  numerical 
stability  can  be  maintained.  The  Spitzer  like  term  is  treated  implicitly.  The  implicit  treatment 
requires  that  the  nonlinear  nature  of  the  diffusion  equation  be  properly  treated,  as  discussed  in 
Sections  II  and  IV.  We  find  that  the  equivalent  diffusion  solution  underestimates  the  effect  of 
preheat  as  compared  to  the  more  accurate  nonlocal  source  solution.  Furthermore,  preheat  as 
calculated  by  equivalent  diffusion  solution  is  susceptible  to  be  further  underestimated  if  there 
is  a  thermal  barrier,  namely  a  region  where  the  actual  energy  flux  goes  up  the  temperature 
gradient.  As  long  as  the  time  step  obeys  the  stability  condition,  the  nonlocal  source  solution 
is  not  held  up  at  all  by  the  thermal  barrier,  and  the  solution  can  be  stably  calculated  even  if  the 
flux  is  up  the  gradient.  However  if  the  nonlocal  source  algorithm  requires  too  small  a  time 
step,  the  equivalent  diffusion  model  can  in  most  cases  be  a  reasonable  substitute. 
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Appendix: 


Here  we  investigate  further  the  physical  reason  for,  and  the  consequences  of  choosing  Scr/T  > 
5.  A  plot  of  W(u)  is  shown  in  Fig  (Al).  Notice  that  for  u  <  3.5,  W(u)  <  0.  This  is  a 
consequence  of  the  fact  that  as  energetic  particles  carry  heat  and  current  forward,  the  less 
energetic  particles  must  carry  current  backwards  so  that  the  total  current  is  zero. 

Hence  if  Scr/T  <  3.5,  the  local  diffusion  coefficient  is  negative.  While  Eq.  (10)  indicates  that 
it  is  possible  that  the  algorithm  be  stable  in  this  case,  we  find  that  in  practice,  this  is  rarely  the 
case.  We  have  solved  the  50keV  slug  problem  using  the  nonlocal  source  solution  and  taking 
Scr/T  =  3.  It  turns  out  that  there  is  a  large  enough  region  in  space  where  the  local  diffusion 
coefficient  is  sufficiently  negative  that  the  code  crashes.  Thus  in  practice,  using  the  nonlocal 
source  solution,  to  maintain  numerical  stability,  we  cannot  let  Scr/T  go  below  3.5. 

The  goal  then  is  to  take  the  minimum  value  of  Scr/T  as  small  as  possible,  so  that  as  many 
particles  as  possible  can  be  treated  nonlocally.  However  we  must  maintain  numerical 
stability  also.  To  us,  a  minimum  value  of  5  seems  a  reasonable  choice.  The  result  is  not  very 
sensitive  to  this  choice  as  long  as  Scr/T  remains  small  and  numerical  stability  is  maintained. 

In  Fig.  (A2)  are  shown  solutions  to  the  50  keV  slug  problem  at  time  t  =  3ns  for  several  cases, 
Krook  model  nonlocal  source  solution  with  Scr/T  =4  and  5,  Spitzer  conductivity,  and  flux 
limited  diffusion  with  a  flux  limit  of  0.06.  Notice  that  even  for  Scr/T  as  small  as  4,  stability  is 
maintained  and  the  solution  is  not  very  different  from  that  with  Scr>5.  These  two  solutions 
are  much  closer  to  one  another  than  they  are  to  either  classical  or  flux  limited  transport. 
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Figure  captions: 

Figure  1 :  A  plot  of  the  stable  region  as  a  function  of  P  (=AAtex/W)  where  At  ex  is  the  explicit 
stable  time  step.  Notice  that  if  P<1  the  numerical  algorithm  is  always  stable. 

Figure  2:  Solution  for  the  temperature  as  a  function  of  time  and  space  using  only  classical 
conductivity.  The  heat  front  advances  quickly,  but  there  is  no  preheating.  The  time  between 
the  curves  is  600ps. 

Figure  3:  a)  As  in  Fig.  2,  but  using  the  equivalent  diffusion  solution  for  the  Krook  model. 
Notice  that  the  main  front  now  advances  more  slowly  and  there  is  some  preheat,  b)  As  in  Fig 
2,  but  using  the  nonlocal  source  solution.  There  is  now  noticeably  more  preheat,  c)  As  in 
Figs  (3  a  and  b)  but  where  the  nonlocal  contribution  to  the  flux  is  artificially  suppressed. 

Figure  4:  a)  The  Krook  model  using  an  equivalent  diffusion  model  where  the  low 
temperature  part  of  the  profile  has  an  initial  temperature  dip.  Notice  that  the  heat  front  has 
difficulty  getting  past  the  temperature  dip.  b)  Using  the  equivalent  source  solution.  Notice 
that  the  temperature  profile  easily  glides  over  the  temperature  dip. 

Figure  5:  A  blown  up  version  of  Fig.  4  showing  a)  the  difficulty  the  front  has  advancing  past 
the  temperature  dip  using  the  equivalent  diffusion  solution,  and  b)  the  ease  with  which  the 
temperature  front  sweeps  over  the  dip  using  the  nonlocal  source  solution.  The  fact  that  the 
flux  is  up  the  gradient  does  not  give  rise  to  any  numerical  difficulty. 

Figure  6:  A  blow  up  like  Fig  5.  for  the  equivalent  diffusion  model  with  and  without  an  initial 
temperature  dip.  The  graphs  are  at  times  0.29,  1.49,  2.69,  and  3.88  ns.  Notice  that  there  is  no 
preheat  in  the  positive  gradient  region. 

Figure  7:  A  plot  of  the  position  of  the  advancing  temperature  front  as  a  function  of  time  for 
two  different  grids,  128  and  256  cells. 

Figure  8:  A  plot  of  the  temperature  profde  using  the  equivalent  diffusion  solution  (solid)  and 
nonlocal  source  solution  (dashed)  for  the  plasma  studied  in  Refs.  (2  and  10).  The  single  solid 
curve  is  the  initial  temperature  profile. 

Figure  A1 :  A  plot  of  the  function  W(u)  as  a  function  of  u.  Notice  that  for  u  <  3.5,  W(u)  <  0. 

Figure  A2:  A  plot  of  the  temperature  as  a  function  of  position  at  a  time  of  3  ns  for  the  50  keV 
slug  problem  discussed  in  Sec  IV.  We  examine  the  effect  of  the  choice  of  Umin  (=Scr/T). 
Notice  that  while  the  solutions  for  Umm  =  4  and  5  do  not  exactly  overlap,  they  are  much  closer 
to  one  another  than  they  are  to  either  classical  or  flux  limited  transport. 
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